In [1]:
import numpy as np
In [2]:
def gauss_elim_nopiv(n, matrix, b):
origin_m = np.copy(matrix)
origin_b = np.copy(b)
aux_m = np.array([[0.]*n]*n)
for k in range(0, n):
# print(origin_m, origin_b)
for i in range(k+1, n):
aux_m[i, k] = -origin_m[i, k] / origin_m[k, k]
origin_m[i, k] = 0
origin_b[i] = origin_b[i] + aux_m[i, k] * origin_b[k]
for j in range(k + 1, n):
origin_m[i, j] = origin_m[i, j] + aux_m[i, k] * origin_m[k, j]
determ = origin_m[0, 0]
for d in range(1, n):
determ *= origin_m[d, d]
return origin_m, aux_m, origin_b, determ
In [3]:
def backward_sub(n, U, c):
x = np.zeros_like(c)
x[n-1] = c[n-1]/U[n-1, n-1]
for k in range(n-2, -1, -1):
s = 0
for j in range(k+1, n):
# print(k, j)
s += U[k, j] * x[j]
# print(k, j, U[k, j], x[j], s)
x[k] = (1/U[k, k]) * (c[k] - s)
# print(x[k], x)
return x
In [4]:
def forward_sub(n, L, b):
y = np.zeros_like(b)
y[0] = b[1] / L[1, 1]
for k in range(1, n):
soma = 0
for j in range(0, k - 1): soma += L[k, j] * y[j]
y[k] = (b[k] - soma) / L[k, k]
return y
In [5]:
def gauss_elim_partial(n, matrix, b):
origin_m = np.copy(matrix)
origin_b = np.copy(b)
for k in range(0, n - 1):
maxi = 0
km = k
for k1 in range(k, n):
if abs(origin_m[k1, k]) > maxi:
maxi = abs(origin_m[k1, k])
km = k1
if km != k:
origin_m[k, k:n-1], origin_m[km, k:n-1] = origin_m[km, k:n-1], origin_m[k, k:n-1]
origin_b[k], origin_b[km] = origin_b[km], origin_b[k]
U_A, aux_A, c_A, det_A = gauss_elim_nopiv(n, A, b)
return U_A, aux_A, c_A, det_A
In [6]:
def solve_nopiv(matrix, b):
n = len(matrix)
u, aux, c, det = gauss_elim_nopiv(n, matrix, b)
x = backward_sub(n, u, c)
return x
def solve_partpiv(matrix, b):
n = len(matrix)
u, aux, c, det = gauss_elim_partial(n, matrix, b)
x = backward_sub(n, u, c)
return x
In [7]:
A = np.array([[1, 2, 1], [2, 3, 3], [-1, -3, 1]])
b = np.array([1, 3, 2])
A, b
Out[7]:
In [8]:
gauss_elim_partial(3, A, b)
Out[8]:
In [9]:
gauss_elim_nopiv(3, A, b)
Out[9]:
In [10]:
U_A, aux_A, c_A, det_A = gauss_elim_nopiv(3, A, b)
print(U_A, aux_A, c_A, det_A)
In [11]:
x_A = backward_sub(3, U_A, c_A)
print(x_A)
In [12]:
# y_A = forward_sub(3, aux_A, c_A)
# print(y_A)
In [13]:
def hilbert_gen(n):
H = np.array([[1.]*n]*n)
H_b = []
for b in range(1, n+1): H_b.append(float(b))
for i in range(0, n):
for j in range(0, n):
H[i, j] = (1/((i + 1) + (j + 1) - 1))
return H, np.array(H_b)
In [14]:
Hil_3, Hil_3_b = hilbert_gen(3)
print(Hil_3, Hil_3_b)
print(np.linalg.solve(Hil_3, Hil_3_b))
print(np.linalg.det(Hil_3))
In [15]:
H3_S, H3_A, H3_B, H3_D = gauss_elim_nopiv(3, Hil_3, Hil_3_b)
backward_sub(3, H3_S, H3_B)
Out[15]:
In [16]:
print(gauss_elim_nopiv(3, Hil_3, Hil_3_b))
In [17]:
Hil_5, Hil_5_b = hilbert_gen(5)
gauss_elim_nopiv(5, Hil_5, Hil_5_b)
Out[17]:
In [18]:
Hil_10, Hil_10_b = hilbert_gen(10)
print(Hil_10, Hil_10_b)
print(np.linalg.solve(Hil_10, Hil_10_b))
print(np.linalg.det(Hil_10))
In [19]:
gauss_elim_nopiv(10, Hil_10, Hil_10_b)
Out[19]:
O determinante deu 2.164480500135833e-53, que é parecido com o obtido pelo cálculo do determinante utilizando o NumPy
In [20]:
A = np.array([[1, 2, 1], [2, 3, 3], [-1, -3, 1]])
b = np.array([1, 3, 2])
A, b
Out[20]:
In [21]:
def matrix_inverse(matrix):
n = len(matrix)
inv_matrix = [[] for _ in range(n)]
for i in range(0, n):
e = [0, 0, 0]
e[i] = 1
U, aux, c, det = gauss_elim_nopiv(n, matrix, e)
x = backward_sub(n, U, c)
inv_matrix[i] = x
inv_matrix = np.array(inv_matrix)
return inv_matrix.transpose()
In [22]:
A_inv = matrix_inverse(A)
print(A_inv)
print(A.dot(A_inv))
Utilizando o NumPy determina-se a inversa de A como:
In [23]:
np.linalg.inv(A)
Out[23]:
Exercício 4 - Resolver um sistema para entender a necessidade do pivoting
In [24]:
C = [[-10e-20, 1], [2, 1]]
cs = [1, 0]
In [25]:
solve_partpiv(C, cs)
Out[25]:
In [26]:
x1 = -0.4999975
x2 = 0.999995
In [27]:
C = [[1, 2], [3, 4]]
D = [[1, 1], [3, 2]]